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1. 


Introduction 


A computational  technique  is  described  in  the  present  paper, 
to  evaluate  steady,  inviscid,  compressible  flows  past  certain 
vehicles  flying  at  supersonic  speed.  The  object  of  the 
analysis  is  twofold: 

1)  to  provide  a tool  for  an  accurate  evaluation  of 
inviscid  blunt  body  flows,  even  in  the  presence  of  imbedded 
shocks,  and  of  flows  about  afterbodies  whose  cross-sections 
are  not  far  from  circles  but  whose  longitudinal  sections 

are  substantially  different  from  sections  of  cones  or  cylinders 
( a typical  situation  occurring  in  vehicles  originally  designed 
as  cones,  bicones,  cones -cylinders , as  a consequence  of 
ablation) ; and 

2)  to  provide  a reliable  basic  computational  tool  for 
a further  extension  to  the  analysis  of  viscous  flows. 

Since  departure  from  axially  symmetric  geometries  is 
assumed  to  be  not  too  strong  (as  opposed  to  geometries  of 
aircraft,  see  Ref.  1),  we  will  stipulate  to  build  our  analysis 
around  a basic  cylindrical  frame  of  reference,  whose  axis 
(to  be  called  x)  lies  in  the  general  direction  of  the  length 
of  the  vehicle.  Meridional  planes,  defined  by  an  angle  cp, 
contain  identical  flows  if  the  geometry  is  axially  symmetric 
and  the  angle  of  attack  is  zero.  Three-dimensional  flows 
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will  be  analyzed  by  adding  three-dimensional  effects  (along 
the  cp-coordinate)  to  a basic  axisymmetric  analysis.  The 


latter  will  thus  be  considered  our  initial  goal,  for  the 
sake  of  simplicity,  and  will  be  described  in  the  present 
paper.  We  will  also  limit  the  analysis  to  flows  without 
imbedded  shocks.  Extensions  of  the  technique  to  flow  with 
shocks  and  three-dimensional  flows  will  be  discussed  elsewhere. 


2 . Generalities  and  Nomenclature 

Let  x be  the  axis  of  symmetry  of  the  vehicle,  and  y 
the  radial  coordinate  in  any  meridional  plane  (x  and  y here 
stand  for  the  more  familiar  z and  r of  a cylindrical  coord- 

A A 

inate  system) . Let  I and  J be  the  unit  vectors  m the  x- 

and  y-direction , respectively . As  usual  (Ref.  1),  we  will 

choose  P=ln(p/pa_))  and  S as  thermodynamical  unknowns  (where 

S is  the  entropy  measured  from  its  value  at  infinity,  divided 

by  c ) and 
v 

V = UI  + VJ  (1) 


to  represent  the  unknown  velocity  vector  which,  as  well  as  its  modu- 

l. 


lus  q,  will 
addition, 
and  of  the 


be  considered  in  multiples 
we  will  make  use  of  Y,  the 

variable,  nondimens ional 
P/P, 


,/  = 


CD 


P/P. 


CO 


of  (p  /:  ) In 

oo  oo 

ratio  of  specific  heats, 
temperature , 

(2) 


A 


~> 


Lenqths  will  be  scaled  to  a reference  value,  x typical  of 

ret  c 

the  vehicle. 


Two  problems  will  be  examined: 

1)  the  evaluation  of  the  steady,  transonic  flow  past 
the  blunt  nose  of  the  vehicle,  for  which  a time-dependent 
computational  technique  will  be  used,  and 

2)  the  evaluation  of  the  steady,  supersonic  flow  about 
the  afterbody,  for  which  a marching  computational  technique 
will  be  used. 

When  time-dependent  calculations  are  performed,  all  times 

p 

will  be  scaled  to  a reference  time,  t =x  ,/(p  /o  ) 

ref  ref  oo  co 

The  equations  of  motion  are,  then,  conveniently  written 
in  the  form: 

- v 

t y I 

Vfc+i2Vq2  -VxtxV+7  P = 0 | (3) 

st+v-vs  = 0 ' 

where  j-1  for  axisymmetric  flows.  The  same  system  may  be 
used  for  two-dimensional  problems  by  simply  letting  j=0.  It 
is  also  clear  that  steady  problems  are  defined  by  (3)  after 
dropping  all  time  derivatives. 
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3 . Computational  Frame  and  Basic  Equations  of  Motion 


Most  of  the  calculations  of  blunt  body  and  afterbody 
flows  have  been  performed,  in  the  last  decade,  using  a 
cylindrical  frame  such  as  the  one  described  above,  or  a spher- 
ical frame,  and  in  many  cases  the  results  have  been  excellent. 
Calculations  about  bodies  with  concavities,  such  as  ablated 
cones  or  bodies  with  flares,  may  not  be  equally  successful, 
primarily  because  the  computational  grid  cannot  be  properly 
adjusted  to  the  body  geometry.  It  is  well-known,  indeed,  that 
one  of  the  grid  lines  should  intersect  the  body  under  an 
angle  as  close  as  possible  to  a right  angle,  in  order  to  im- 
prove accuracy. 

Therefore,  following  the  same  line  of  thought  as  exposed 
in  Ref.  1,  but  in  meridional  planes  rather  than  in  cross- 
sectional  planes,  we  will  recast  the  equations  of  motion  in 
another,  curvilinear,  frame  defined  by  two  variables,  l and 
T|,  such  that  the  image  of  the  body  in  the  ( % , 1)  -plane  be  as 
close  as  possible  to  an  T|-constant  line.  To  this  effect,  we 
will  use  a conformal  mapping  to  represent  the  correspondence 
between  the  complex  z-plane,  where 

z = x+iy  (4) 

and  the  complex  1 -plane,  where 

C = i+i1!  (5) 
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r 


Details  on  the  nature  of  the  conformal  mapping  will  be 
given  later  on  (in  Section  14) . We  can  provide,  however,  a 
large  number  of  general  formulas,  valid  regardless  of  the 
particular  mapping  being  used. 

A A 

Let  i and  j be  the  unit  vectors  (in  the  z-plane)  parallel 
to  the  %constant  line  and  to  the  ?=constant  line,  respectively. 
Let  us  also  define  new  velocity  components,  u and  v,  by 


V 


A A 

ui+vj 


(6) 


Let  us  also  define 


g = 


cT 

dz 


Ge 


ix 


1_  d log  g 
g dz 


cpi  +icp2 


and 


It  is  convenient 


C+ i£  = 

to  note 


G -iuu 
— = e 

g 

that 


(7) 

(8) 

(9) 


~ =GC  , % =G>  , 
-x  -y 

T|  =-GS  , D =G C 
x y 

(10) 

C/G  , x_=-£/G  , 

y - =£/g  . yr=C/G 

(ii) 

lZ+^=Gs  , 

xl+x2  =1/GS 

(12) 

x y 

i V 

:Gcpi  / Gp=-Gcp2x  , 

'-'r=CPp  > 

(13) 

r * x c * 

= C , 1 • J-A-  , j • 

A A A p 

i =-£  , j-J=C 

(14) 

U=uC-v£  , v=u £+vC  , u=UC+v£  , v=-U£  +VC  (15) 
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The  basic  vector  operators  in  (3)  can  be  expressed  as  follows: 


Vp=G (P_i+P?3) 


■ xV=G2  [(r-)  _-£)  T]h 
G : G 'I 

— _ — - r , v,  ,u.  -i  , A A, 

Vx  -'xV=G ' - (-)  y-(~)  (V1-U3) 


(16) 


so  that  (3)  become: 

P_+GuP_+GvP  ^4*  yGu_+ YGv„  — Y ( uG _ + vG -, ) + i \ = 0 

T - T)  e M c,  U y 

u,+Guu,+Gvu^+vsG.-uvG.+  <JGP_  = 0 
T ? M i M 

v +Guv_4-GvvT1-uvG_+usGT1+  (7GP-  = 0 

T ; ; I,  1| 

ST+G  (uS_.+vS?i)  = 0 

where  " has  been  used  instead  of  t,  although  derivatives 
with  respect  to  t and  with  respect  to  t are  identical,  since 
the  mapping  does  not  depend  on  time. 

We  will  now  make  use  of  (13)  and  introduce  the  notations: 

V 

D=vo1+ucpE  , E=-ucp! +v¥;;  + j — (18) 


to  rewrite  (17)  in  the  simpler  form: 


P_+G  (uP.+vP^+Vu.  + Wp+vE)  = 

u +G(uu.+vu„+vD+  2P.)  = 0 

T =,  M • 

v +G  (uv_  + vv„-uD+  2P_)  = 0 
•I  M 

ST+G(uS_  + vSj  = 0 

1 " M 


(19) 
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t . Coordinate  Stretching  and  Final  Equations  of  Motion 

In  viscous  flow  problems,  a strong  concentration  of  grid 
lines  near  the  body  is  necessary,  to  describe  the  boundary 
layer  with  sufficient  resolution,  and  it  can  be  achieved  by 
stretching  the  T( -coordinate . We  will  define  a stretching 
function  depending  on  one  parameter,  a ; values  of  a less  than 
1 produce  almost  undetectable  stretchings,  whereas  the  concen- 
tration of  grid  lines  near  the  body  can  be  made  as  great  as 
necessary  by  using  sufficiently  large  values  of  a. 

In  general,  let 

Tl  = b ( f ) (20) 

and 


T\  = c ("  , t) 


(21) 


be  the  equations  of  the  images  of  the  body  and  of  the  shock 
in  the  '-plane.  We  will  define  a new  coordinate,  X,  as 
follows : 


X=l+  — argtanh  , 'n=c+$tanh[u  (X-l)  ] (22) 

a & 


where 


& = 


c -b 


tanh 


(2  3) 


Therefore, 


X_= 


Tl  a[&2  — ( li  — c ) ] 


(24) 


1 


7 


,-c 


-(c_-b  J+c_]x, 


X?  '&tanh 


X =-[— — +l]c_X 
t Stanh  a 


(25) 

(26) 


Letting  Y=- , T=t , 


A = X +GuX_+GvX, 


(27) 


the  equations  of  motion  (19)  become 


pt+apx+g(uPy+Vuy+yxsux+vxtivx+ye)  = ° 

u^+Au^+G  (uu^,+  vD+  ??y+  7X_P  . ) = 0 
v +Av  +G (uv  -uD+  ?X_P  ) = 0 

1 A X 1 1 A 


ST+ASx+GuSY  = 0 


Note  that,  at  the  body,  71=b,  X=0,  X„==cosh' , 
L.  = -b_XT,,  X.=0,  and,  at  the  shock,  7=c,  X=l,  X_=1/o.jS, 


(28) 


X_  = 


-c.Xr  X = -c  X_ . 

§ U T t u 


In  general, 

?x=o,  ry.l,  ;T=o,  V1/Xr  V’W  V'Vb 


(29) 


5 . Characteristic  Equation  for  the  Time-Dependent  Problem 
To  evaluate  points  on  the  body  surface  or  on  the  bow 
shock,  a characteristic  equations  in  X and  T,  obtained  from 
the  first  three  equations  (28),  is  needed.  With 

Ri =uPy+Yuy+\E 

P.s=uuy+vD+  7 P Y } (30) 

Rt  =uvy-uD 


1 


a linear  combination  of  the  first  three  equations  (28)  is 
written  in  the  form: 

3 

ui  (PT+\px)+u2  (ut+1ux)  + u3  (vt+:'-vx)+G  ^ 4iRi=  0 (31) 

1 

whose  X turns  out  to  be  defined  by 

\=A±aG'/x3+X3  (32) 

and,  consequently. 


Ui  = A-X , Us  =-YGX_ , u3=-VGXt. 


(33) 


6 . Body  Calculation  for  the  Time -Dependent  Problem 

At  body  points,  X is  defined  by  the  lower  sign  in  (32) 
and  A=0.  The  boundary  condition  yields: 

V • N = uNx+vNs  =0  (34) 


AAA 

where  N=Nii+N2j  is  the  unit  vector  normal  to  the  body. 

A A 

* _ -dyi+dxJ 
(dxc +dy3 ) ^ 


Since 

(35) 


(11),  (12),  (14)  and  (20)  can  be  applied,  to  obtain 

-dyI+dxJ=-  (y_dP+y^dTl)  1+  (x^dc+x-dfl)  J 

= J C-(id?+CdTl)  1+  (Cd?-£dT|)  J] 

1 r A-, 

= - [ -d  'ii+d?  j ] 

N2=  ^ , Ni  = -b_N2  , v=vl+b3  (36) 

■ 
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It  follows  that 


-ub_+v  = 0 


-u  b.+v  = 0 
T ? T 


By  replacing  (33)  into  (31)  and  taking  (38)  into  account. 


one  obtains: 


-X(PT+\Px)  + :vVXT1(bEux-vx)  + h„iRi=0 


Therefore,  the  pressure  is  defined  by: 


♦ t "iRi 


The  entropy  is  obviously  constant  at  the  body,  if  initially 
assumed  equal  to  the  stagnation  entropy.  The  velocity  can 
be  computed  as  follows:  Let 


v = u-rvb_=uv< 


Multiply  the  third  of  (28)  by  b_  and  add  to  the  second. 


At  the  body. 


which,  since 


UT+VTbr  G'U(UY+VYbr.)+  7PY:  = ° 


v = u +v  b 
T T T - 


VY=uY+VYhr+vb?- 


can  be  written  as 


vyf-Giu (v  -vb._) + 7Py]  = 0 
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(42) 


(43) 


yields 


-A..+ C’_  y'  + Cy  x.-  ~ ( C-+X-V  ’ +£y  x_)  2/ (C+S-y' ) 

S -3  3 CJ  “ -5  C 


dr 


so  that,  by  taking  (9),  (11)  and  (13)  into  account, 

( 1+ y 7 2 ) ^.?  +Cy " / g 


d°  b 
d--? 


(C+^y') 


(44) 


7 . Bow  Shock  Calculation  for  the  Time-Dependent  Problem 

To  compute  the  bow  shock  in  the  time-dependent  problem, 
we  define  the  unit  vector  normal  to  the  shock,  N,  as 


AAA 
N = N;  i + Nr  j 


(45) 


By  arguments  similar  to  the  ones  leading  to  (36) , we  obtain: 

Ns  = l/v  , Ni=-c_Nr  , Vssv'l+c!.  (46) 

Note,  however,  that  c depends  on  t as  well  as  on  f.  Therefore, 
vT=c?c=T/v  = cycYT/v  (47) 
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N2t=-vt/V  , N.  T— cyt/v+CyVt/v‘ 


The  free-streara  velocity,  for  an  axisymmetric  flow,  must  be 
parallel  to  the  x-axis.  Consequently, 


u =U  r 

oo  oo c 


U = -V  - 
coT  ao  T 


v = -U  £ 

OO  CO 


v =u  r 
oo  T oo  T 


If  Q is  a point  on  the  shock,  the  velocity  of  the  shock,  W, 
is  given  by 

dQ  A . a a a A 

W = ~~  -N  = (x  I-N+ (y  ) j-n  ( 

dt  Q t Q t 


(XQ,t+i(yQ,t  ‘ <Vt*  g<:Q't 


Now,  let  Q be  a point  which  remains  on  the  same  inconstant 
line,  as  the  shock  moves.  Therefore, 


(Vt+i,Vt=  g Ct 


From  (52)  and  (50) , using  (45)  and  (14) , we  obtain 


W = — Nrc 
G T 


Let  u be  the  velocity  component  normal  to  the  shock,  and 

u , = u -W 
rel 

Consequently , 


u = u N;+(v  - — c ) N- 

oorel  oo  oo  GT 


and 


(u  , ) =-v  Ni  a +u  Ni  T+ (u  -u  + TY*  G_c  - ^ cmm)  Na+ (v  - ~ c 
oorel  T oo  T oo  1 cd  T G TTG  TT  co  G 


Note  that  * =0  but 


X 

t 

uu  = -cp,  — = Cpi  c 
T 1 X?  * T 


at  the  shock.  Similarly, 


(56) 


(57) 


GT=-cp2cTG 


(58) 


We  will  write  (56)  in  the  form: 


by  letting 


(u  _ ) =C'  c +Cj 
core!  T * TT 


Ci  =-N2/G 


C2  = (u  N2 -v  Ni ) w +u  NiT+(v  - — c_)  N2rn+cp2c*Lc 
co  oo  T co  T on  r,  T'  ‘T  t 


oo  G T 


(59) 

(60) 
(61) 


We  write  now  the  Rankine-Hugoniot  conditions  in  the  form: 

2 , Y-l, 


P = In 


, + tn(ur 
Y+l  oorel  2 


(62) 


Y-l 


2 Y 


rel  Y+l  oorel  Y+l  u 

oo  rel 

By  differentiation  with  respect  to  T,  and  by  letting 

2u 


(63) 


Di  = 


oo  rel 


u 


Y-l 


, Dr  = 


Y-l  _ __2Y  1_ 

Y + l Y + l uT 


oo  rel 


oo  rel 
(64) 
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T)  N2T 


J 


we  obtain  from  (6  2)  and  (6  3)  : 


with 


P = C3  c +C4 
T TT 


(u  ,)  = C_-c  +C... 

rel  T TT 


(65) 

(66) 


with 


Cl  • = C9  N; 


Cx3=C9Na 


C-  2 = -V  i+Ci  ;N-  + (u-u  ) Ni  rr 
oo  T ' OD  • T 


C-  i=u  1 +C-  (u-u  ) N 

oo  T oo  1 


(75) 


We  use  now  the  characteristic  equation  (31)  with  V defined 
by  (32)  with  the  upper  sign,  and  consequently 


lli  — A—  X.  , b*2~7GX_c.  , us  — — YGX  — 


(76) 


By  replacing  (74),  (75)  and  (65)  into  it.  and  letting 

D?  =1/  (di  Ca  + u^Ci  i + i_i3  Ci  3 ) 

Dg  =Ux  C-1  + U.2  Cl  C + U3  Cl  •'» 


Dg  =GHu . R . 

1 1 


n=Mux Px+uaux+a3 v ) 


the  following  equation  is  obtained,  to  define  c : 


(77) 


ctt=-d7  (D8+D3  + n) 


(78) 


8 . Outline  of  the  Computational  Program  for  the  Time-Dependent 
Pjroblem 

The  steady  flow  about  the  blunt  nose  of  the  vehicle  is 
determined  as  the  asymptotic  solution  of  an  unsteady  flow 
problem.  A computational  region  is  chosen,  subject  to  the 
only  condition  that  the  sonic  line,  as  anticipated  in  the 
steady  state,  must  lie  entirely  inside  it.  A guess  is  made 
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for  the  initial  standoff  distance  of  the  Low  shock,  d , as 

o 

a multiple  of  the  nose  radius,  as  follows: 

d = 2.7687/ (M3  -l)+.36  for  two-dimensional  (79) 
o co 

flows 


d = 0.6137/(M2  -1)+.13  for  axisymmetric 

o CD  .... 

flows 


(80) 


After  evaluating  the  image  of  the  normal  shock  point 
in  the  '-plane,  the  other  shock  points  are  assumed  to  have 
the  same  ^-coordinate  as  the  former.  A computational  grid 
can  be  defined  now,  and  the  shape  of  the  shock  in  the  physical 
plane  can  be  found.  We  assume  that  the  shock  is  initially  at 
rest;  therefore,  cT=0,  17=0.  The  normal  to  the  shock  at  each 
point  is  defined  by  (45)  and  (46) , the  free  stream  velocity 
components  by  the  first  set  in  (49)  and  (55)  can  be  used, 
with  cT=0.  Pressure  and  velocity  behind  the  initial  shock 
are  then  defined  by  (62) , (63)  and  (73) , whereas  the  entropy 

follows  from 

S = P-Y'.n(u  /u)  (81) 

oo  rel 

The  entropy  on  the  body  surface  is  assumed  equal  to  the 
stagnation  entropy  and  initial  guesses  for  P and  the  modulus 
of  the  velocity  on  the  body  are  made,  assuming  a Mach  number 
distribution  along  the  body  and  using  well-known  expressions 
for  steady,  isentropic  flow.  The  two  velocity  components 
along  the  body  are  obtained  by  applying  the  kinematic  boundary 
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condition,  that  is,  vanishing  of  the  normal  velocity  component. 


1 


Values  of  P,  S,  u and  v at  all  other  grid  points  are  initially 
prescribed  by  linearly  interpolating  between  body  and  shock. 

Each  computational  step  consists  of  a predictor  level  and 
a corrector  level. 

Predictor  level. 

Given  original  values  of  P,  S,  u,  v,  7,  c,  c , z,  ',  g 
and  ip, 

At  all  values  of  ?,  compute  c^,  c . 

At  all  values  of  ? and  "H, 

compute  X-  and  Y-derivat ives  of  P,  S,  u,  v (see  next  section) 

compute  X_,  X , X 
-'It 

compute  G,  A,  D,  E. 

For  all  points  except  body  and  shock  points,  determine  P , 

S , u and  v from  (28) . 

T T T 

For  body  points,  compute  R.  , R?  , Ra  from  (30)  , \ from  (32)  , 
and  u:  , u:  , it-  from  (33).  Then,  use  (39)  to  obtain  P , (41) 
to  obtain  and  set  S^=0. 

For  shock  points,  compute  R-.  , R-  , Rj  from  (30)  , evaluate 
Ni  and  Nr  from  (46)  and  compute  all  coefficients  (C-.  through 
Ci -i ) . Compute  D~  , D3  , Dg  and  and  determine  c from  (78)  . 

Update  P,  S,  u,  v at  all  interior  points,  P and  v at  all 
body  points,  as  well  as  cT  and  c,  using  the  following  rule 
(where  ‘ is  an  arbitrary  function: 
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(T-f-iT) 


(32) 


= ?(T)  + ?t:-T 

From  the  new  geometry  of  the  shock  (determined  by  the 

updated  values  of  c)  determine  the  basic  parameters  for  a new 

mapping.  From  the  updated  values  of  c,  determine  new  values 

of  c . and  new  values  of  Ni  , N2 . Then,  new  values  of  u 

Y ' co  rel 

are  determined  and  (62) , (63) , (73) , (81)  used  to  evaluate 

P,  u,  v,  S behind  the  shock.  Values  of  .7  are  now  computed 
throughout  by  using: 

7 = exp(~ ^ P+  ~ S)  (83) 

Symmetry  conditions  are  imposed  where  necessary.  The  initial 
values  of  P,  S,  u,  v are  temporarily  saved. 

A new  grid  (consistent  with  the  new  shock  location)  is 
computed,  including  the  values  of  g and  cp. 

Corrector  level. 

The  computation  is  restarted  as  at  the  beginning  of  the 
predictor  stage,  using  all  new  values.  The  updating  is  now 
performed  using  the  following  rule  in  lieu  of  (82)  : 

£(T+iT)=  (T)+? (T+1T)+?T1T]  (84) 

where  ? (T+1T)  is  the  left-hand  side  of  (82)  and  ? is  the 
T-derivative  computed  in  the  corrector  stage.  The  computa- 
tional grid  must  be  re-evaluated  since  the  bow  shock  location 
may  have  changed  slightly.  The  values  at  the  shock  them- 
selves have  to  be  recomputed  as  at  the  end  of  the  predictor 
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level.  Symmetry  conditions  are  imposed  again.  Finally,  the 
updated  values  of  P,  S,  u,  v,  7,  c and  cT  are  stored  as  initial 
values  for  a new  step. 

9 . Discretization  of  X-  and  Y -Derivatives 

In  principle,  the  MacCormack  scheme  (Ref.  2)  for  inte- 
grating the  equations  of  motion  at  all  interior  points  is 
adopted.  This  policy  is  reflected  by  the  use  of  (32)  and  (84) 
at  the  predictor  and  corrector  level,  respectively.  In  addi- 
tion, whenever  possible,  the  X-  and  Y-derivatives  should  be 
approximated  by  2-point  differences  taken,  for  example,  forwards 
at  the  predictor  level  and  backwards  at  the  corrector  level. 
There  are  several  exceptions  to  the  above  rule: 

1)  At  body  boundary  points,  where  X=0,  backward  approx- 
imations cannot  be  taken  and  are  replaced  by  the  approximation: 

$ *(-25  +3?!  -i2)/lX  (85) 

X o 

which  maintains  a second  order  accuracy  if  the  equations  are 
linear  throughout  the  integration  step  (Ref.  3) . Here, 

;.=$ (iAX, Y) . 

2)  Similarly,  at  shock  points,  where  X=l,  forward  approx- 
imations cannot  be  taken  and  are  replaced  by  the  approximation: 

« «(2i  -3?;  + iE)/.'X  (86) 

X o 

where  1 =?  (1-i; X,  AY)  . 
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3)  The  last  of  (28)  simply  expresses  the  fact  that 
entropy  is  transported  unaltered  by  the  moving  particles. 

That  means  that  no  entropy  signal  travels  backwards  along  a 
particle  path.  Consequently,  no  approximations  to  X-  or  Y- 
derivatives  are  allowed  which  imply  forward  differences  of 
S along  a particle  path.  Violation  of  this  physical  principle, 
as  a consequence  of  straightforward  application  of  the 
MacCormack  scheme,  affects  the  results  along  the  symmetry 
streamline  in  a blunt  body  problem,  and  neighboring  points 
are  also  affected  consequently , with  an  overall  unnecessary 
lengthening  of  the  computational  time  for  reaching  a steady 
state  (Ref.  4) . It  is  therefore  advisable  to  adopt  the 
following  scheme  systematically: 


Predictor  level:  S » sgnA 

X. 


Sy~  sgnu 


2So o - 3S^o+S2o 
IX 

2 So  O —3 So  1 + S; 2 
A Y 


So  o -Si o 


Corrector  level:  S^x:  sgnA  — 7 ~ 


So  0 - Sq 1 


~ sgnu  — — 


where 

So  0 =S (X, Y)  , Si  0 =S (X-AXsgnA,  Y) 

Soo=S  (X-2  AXsgnA,  Y) 

Sc  1 =S (X, Y -iYsgnu)  , S02 =S (X, Y-2£Ysgnu) 


(87) 


(88) 
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4)  Similarly,  the  terms  v +Guv  in  (41)  express  the 
Lagrangian  derivative  of  v.  Forward  differencing  of  v is 
forbidden  and  v^  should  be  approximated  as  follows: 

Predictor  level:  v^.  ~ ( vo  -v-,  ) / LY  j 


Corrector  level:  v « (2v0  -3v-_  +v£  ) /-Y 


vo=v(o,Y)  , v-i  =v  (o,  Y~LY)  , v2~v (o, Y-21Y) 


5)  For  the  same  reason,  that  is,  for  a better  evalua- 
tion of  Lagrangian  derivatives  in  the  momentum  equations,  the 

values  of  u , v , u , v in  the  second  and  third  of  (28)  should 
X X Y Y 

be  approximated  following  the  same  rule  as  described  by  (87) . 

6)  Y-derivatives  at  the  bow  shock  are  conveniently 
approximated  by  centered  differences,  except  for  c^  at  the 
last  two  points,  where  two-point  backward  differences  are  used 
(Ref.  4) . 


10 . Equations  of  Motion  for  the  Afterbody  Flow 

In  studying  the  steady,  supersonic,  axisymmetric  flow  in 
the  shock  layer,  we  may  start  again  from  (17) , after  dropping 
the  time  derivatives,  and  introducing  an  angular  variable. 


After  some  algebraic  manipulations,  the  following  set  of  equa- 


tions is  obtained: 
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xp_+aP_+Y[a_-(l+~2 ) ~f-  +j  ~~]  = 0 
5 T1  T|  G Gyu 


'*+1Z?t+  ua  + (1+:r  )(Q 


-)  = 0 


(92) 


S.+rS_  = 0 

= 4 


where 


K = 1 


u 


(93) 


These  are  the  only  equations  which  require  integration  since, 
for  a steady,  isenergic  flow, 

2 Y 


* - tT(V7) 


(94) 


where  7 is  the  stagnation  value  of  2,  and,  in  turn, 

Y-l  1 

U=  exp  (~  P+  - S) 


(95) 


If  the  same  stretching  function  as  in  Section  4 is  used,  and 
X is  given  the  same  meaning,  but  T is  now  made  equal  to  - , 
we  may  use  the  following  definitions  in  order  to  maintain  our 
present  notation  as  close  as  possible  to  the  one  used  in 
Ref.  1: 


A = V 7 XT1  ' B = VCXT1 


V 


C=X_-UX.  , D=(l+-B)  CPi -j  , E=(l+0‘  ) (JCPa+cps) 

1 ^ Cjyu 


(95) 

(97) 


Since,  for  any  function,  F: 


F-=VFxx-  ' FrFxxTi 


(98) 
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the  system  (92)  becomes: 


pt+apx+  t <Vx‘D)  = 0 


:T+B'X+  -(CPX-rPT)  'E  = 0 


S+BSV  = 0 
T X 


Eqs.  (29)  must  be  replaced  by 


T] 

X 


't 


= -X  /x„  . 

=.  M 


(99) 


(100) 


11.  Characteristic  Equation  for  the  Steady  Flow  Problem. 

Body  Calculation 

To  evaluate  points  on  the  body  or  on  the  bow  shock,  a 
characteristic  equation  in  X and  T,  obtained  from  the  first 
two  equations  (99) , is  needed,  in  the  form  of  a linear  combination: 


(Ui-3  ^ via)  (PT+iPx)  + u2  (3t+1cx)=  J Du1+E. 


whose  1 is  defined  by 

\=B+  ^i-X-Ca±3]  . e2  = (l+J2)  K-  - 1 

x.u  M a 


and,  consequently , 


4l  — B — X , u.2—  X««  , — - 2 ^2  + £ 

A ' | U A u 


(101) 


(102) 


(103) 


At  the  body,  the  lower  sign  must  be  used;  in  addition,  B=0. 


From  (37) , 


(104) 
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and 


(105) 


Therefore,  the  equation  (101)  defining  the  pressure  at  the 
body  becomes : 


V " XPX+  37(b^+vVE)  ' 7?(J-S)D  (106) 


Obviously,  bT  and  b. _ will  be  evaluated  according  to  (43)  and 
(44)  . 


12.  Bow  Shock  Calculation  for  the  Steady  State  Problem 

Equations  (45)  and  (46)  still  hold,  as  well  as  (49) . 

By  writing  c'  and  c"  in  lieu  of  c_  and  c?t_,  respectively, 
we  have  now: 

v =c'c"A  / NaT  =-  c'c"A3  , Nxt=-  c ' A3  (107) 

The  basic  unknown  is  still  c , that  is,  c" . To  determine 

it,  we  note  first  that 

u = u Nx+v  N2  (108) 

oo  oo  oo 


Then, 


but 


U _ = U „Nx  + V „Ns  +U  N’  rp+  V N2rp 
coT  coT  ooT  co  ' 1 co 

= U (C  Nx  - > N~  - £c"/v3  + £c 'c"  A3 
oo  T T 


_iii'  -in) 

C_+  i&  = -ie  (r_+r  c )=-ie  (v2+c 

li  - 


so  that 
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(109) 


U = (-V  N-  +U  N2  ) (^2  + C ' ) - (u  + c'v  ) c"/v: 

co  T co  ‘ oo  co  oo 


We  will  then  write 


u = C,  c + C; 
oo  T * TT 


(HO) 


with 


C-  =-  (u  +c  ' v ) /v: 
oo  oo 


(111) 


C2=(-v  N-  +u  N2  ) (t's  + c'o-  ) = -C-  v2  (Ps  + c'O:  ) (112) 

oo  oo  * 


Since  (62)  and  (63)  still  hold,  with  u and  u in  lieu  of 

oo 

u , and  u , , respectively,  it  follows  that  Eqs . (64) 

rel  oorel 

through  (68)  are  valid.  Eqs.  (70)  and  (71)  are  also  valid. 

So  are  (72)  and  (73) ; therefore,  (74)  can  still  be  used,  pro- 
vided that  Ci  i , Ci  s , Ci  3 » Ci  4 are  redefined  as  follows: 


Cl  ! =Cg  Ni  -(u-u  )/v» 


Ci3=C9N2-(u-u  )c'/v3 

oo 


C i 2 ~ — V (Cfc  +C  'cpl)+CioNi 

oo 

Ci4=u  (o2  + c vi  ) +Ci  CN: 

co 


(113) 


Consequently, 


' ip  Cl  £ cTT+Ci  3 


(114) 


with 


C 15=  (Cl  3 -~Ci  1 ) /u 
Ci  3 = (Ci  4 -rCi  2 ) /u 


(115) 


We  use  now  the  characteristic  equation  (101)  with  the  upper 
sign  in  the  definition  of  X;  by  letting: 
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D-7=l/(7:C3+C:iUS) 

Da  = Jp  + Ci  a u 

(116) 

D3=u2  (utD/X^-E) 
r =\(?3Px+ua  :y.) 

the  following  eauation  is  obtained,  to  define  c : 

^ - TT 

CTT=_D7  + (117^ 

This  equation  is  formally  the  same  as  (78)  . 

1 3 . Outline  of  the  Computational  Program  for  the  Afterbody 
Flow 

Initial  conditions  for  a marching  integration  technique 
are  the  converged  steady  values  on  the  next  to  the  last 
5=constant  line  of  the  blunt  body  flow  calculation. 

Each  computational  step  consists  of  a predictor  level 
and  a corrector  level. 

Predictor  level 

Given  original  values  of  P,  7,  S,  c,  c , z,  g 
and  cp  (at  a constant  §)  , 

compute  the  X-derivatives  of  P,  c,  S, 
compute  X_  and  X^, 
compute  G,  A,  B,  C,  D,  E. 

For  all  points  except  the  body  point  and  the  shock  point, 
determine  P , c and  S from  (99) . 
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For  body  points,  determine  Pm  from  (106)  and  let  S =0, 

T T 

For  shock  points,  compute  all  the  coefficients  (C-_  through 
CN  o , D~  , Dj  , D9  , tt)  necessary  to  determine  c from  (117). 

Update  P,  S and  " at  all  interior  points,  P and  S at  the 
body  and  c and  c^  using  (82) . 

With  the  new  value  of  T,  that  is,  new  values  of 

b,  b.and  b_  . are  evaluated,  and  an  updated  c at  the  body  is 
taken  equal  to  b_. 

Moreover,  a new  grid  line  is  obtained,  extending  from  b 
to  the  updated  value  of  c.  New  grid  points  are  determined 
and  the  updated  values  of  P,  S and  c are  affixed  to  them. 

At  the  new  shock  point,  the  new  value  of  c is  used  to 
evaluate  a new  normal  and  then  (62),  (63),  (81),  (72)  and  (91) 

determine  P,  S,  and  c behind  the  shock.  Values  of  7 are  ob- 
tained throughout  from  (95)  and  values  of  u from  (94)  and 

u2  = q2/(l+j2  ) (118) 

The  new  values  of  P,  S,  c,  c,  and  c are  temporarily  saved. 

At  the  new  grid  points,  g and  cp  are  computed. 

Corrector  level 

The  computation  is  restarted  as  in  the  predictor  stage, 
using  all  new  values.  The  updating  is  performed  using  (84). 

The  computational  grid  must  be  re-evaluated  after  correcting  c. 
Finally,  the  updated  values  of  P,  S,  r,  c,  and  c are  stored 
as  initial  values  for  a new  step. 


27 


14.  Mapping  Techniaue 


The  preceding  analysis  holds,  regardless  of  the  technique 

used  to  generate  the  computational  grid.  In  the  current 

series  of  tests,  a very  simple  mapping  technique  was  used, 

★ 

suggested  by  D.  Hall  in  a private  communication.  The  tech- 
nique is  easier  to  analyze  and  to  code  than  the  one  originally 
considered,  which  was  a byproduct  of  the  method  I developed 
for  three-dimensional  calculations  of  supersonic  flows  (Refs. 

1 and  5) . The  original  mapping  maintains  the  geometrical 
symmetry  about  the  body  axis,  if  the  body  is  symmetrical, 
whereas  the  present  one  does  not.  If  the  body  is  not  symmet- 
rical, and  for  three-dimensional  problems  in  general,  con- 
servation of  symmetry  trhoughout  the  mapping  seems  not  to  be 
a necessary  requirement.  In  the  current  series  of  tests,  where 
symmetry  is  still  a typical  feature  of  the  flow,  loss  of 
symmetry  in  the  calculation  of  the  blunt  body  flow  seems 
not  to  affect  the  results. 

A number  of  "hinge-points"  are  selected  inside  the  body, 
to  represent  a coarse  skeleton  of  the  body  itself.  Let  J be 
the  total  number  of  hinge-points,  the  first  of  which  lies  on 
the  x-axis,  in  the  vicinity  of  the  body  stagnation  point. 

Let  us  consider  a succession  of  J+l  planes,  called  the  z_.  - 
planes  (j=l,  J+l),  the  first  of  which  is  identical  with 
the  original  physical  plane  (the  z-plane) , so  that 
♦General  Electric  Co.,  Valley  Forge,  Pa. 
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z-  = z 


All  these  planes  are  related  to  each  other  in  a similar  wav. 


that  is,  through  the  general  function: 

5. 


z . = 1+ (z . -h  . 


1 


'j+1 


1 Dl 


(119) 


where  h_.  ^ is  the  image  of  the  j -th  hinge-point  in  the  z_. -plane, 

and  5 . is  defined  as 
D 


6 = 


j 6 -arg (h  . . -h  . . ) 

D+l.D  11 


(120) 


where  h ^ j is  the  image  of  the  (j+1) -th  hinge-point  in  the 
Zj-plane.  Since  all  points  of  interest  (including  the  hinge- 
points)  lie  in  the  upper  half  of  the  z-plane,  the  set  of  hinge- 
points  will  finally  result  aligned  along  the  real  axis  of  the 
Zj+^-plane.  Note  that  the  image  of  the  x-axis  to  the  left  of 
the  first  hinge-point  will  result  into  a portion  of  the  real 
axis  on  the  zT  -plane,  to  the  left  of  the  image  of  the  first 
hinge-point  (hi  j+]_)  . Note  also  that,  with  the  exception  of 
the  first,  straight  segments  in  the  z-plane  between  two 
successive  hinge-points  are  not  mapped  onto  straight  segments 

in  the  z^  . -plane. 

J+1 

For  each  partial  mapping,  let 


9 • = 


dz  z . -1 

JiJL  = £ 


j dz 


j z . -h  . . 
1 DD 


(121) 


A final  mapping  is  then  applied,  of  the  form: 
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■ '2j+rhi'J+1 


(122) 


for  which 


(123) 


In  the  C -plane,  the  image  of  the  body  becomes  a line  defined 
by  an  almost  constant  value  of  starting  at  c=0,  whereas 
the  portion  of  the  x-axis  in  front  of  the  body  becomes  a 
portion  of  the  imaginary  axis.  Prom  (121)  and  (123)  it 


follows  that 


9 = « jii9^ 


(124) 


According  to  the  defnition  (8) , 


1 _ sl^  V J_  ffi 


9j  dz 


i i ri?i  Vr 

C + 9X95  dzj  1.09'- 


(g  =D 

o 


and,  since 


it  follows  that 


6 .-1 


dz . z . -h  . . j 
3 3 3 3 


i -j  J 5.-1 

1 1 f T 

T + ~ '■  T 

c g z . -h  . . 

j = l ^ ” 


(125) 
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15 . Results  of  Test  Calculations 

We  report  here  results  for  the  following  test  calculations: 


Case  No. 

1 

3 

4 

5 

6 

7 

Free  stream  Mach  number 

3 

4 

4 

4 

4 

A 

Flow  ( two-dimensional 
or  axisymmetr ic) 

2-D 

2-D 

A -3 

2-D 

A-S 

AS 

Body  geometry 

parabola 

sphere- 

cylinder 

sphere-cone 

see  Fig. 9 

Location  of  hinge- 
points  in  physical 
plane 

j-0.75 

I 2 

. -0.9 
-0.63  (1-i) 
0. 9i 
3+0. 9i 

-0.9 

-0.63  (1-i) 
2+1. li 
4 

Number  of  intervals  in 
blunt  body  calculation 

5x12 

6x6 

6x6 

6x6 

6x6 

6x7 

Number  of  intervals  in 
afterbody  calculation 

10 

24 

24 

24 

24 

24 

Stretching  is  not  used  in  any  of  the  preceding  cases.  In 
case  No.  8,  the  same  data  are  used  as  in  case  6 but  a strong 
stretching,  with  a=2  and  12  intervals  between  shock  and  body,  is 
used.  Fig.  1 shows  the  grids  as  they  appear  in  the  physical  plane 
in  both  cases,  for  comparison.  Figs.  2 through  8 show  isobars  as 
computed.  Note  that  there  is  no  difference  between  Fig.  6 and 
Fig.  8.  Finally,  Fig.  9 shows  some  lines  from  the  grid  for  case  7, 
to  demonstrate  the  ability  of  the  mapping  to  adjust  to  bodies  of 
unusual  shape. 

All  calculations  were  performed  on  Polytechnic's  IBM  360/65. 
Their  success  is  due  to  the  dedicated  cooperation  of  Miss  Catherine 
Fahy,  without  whose  help  I would  not  have  been  able  to  get  my 
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program  to  operate  on  the  machine's  extremelv  complicated 

system. 
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COMP  ITT  AT  I ON  AL  GRIDS  FOR  RUNS 


FIG.  3.  ISOBARS  FOR  RUN  #3 


FIG „ 5.  ISOBARS  FOR  RUN  #5 
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FIG.  8 - ISOBARS  FOR  RUN  #8 


